---
title: "Replication_Diener_Age_Rep_Styles"
output: html_document
---
# Readme

These are the replication files for the Research Note "Age and Representation Styles: Are Young Candidates more likely to Prioritize their Voters over their Own Opinion or their Party?" by Julius Diener.

The following script contains all code required to reproduce the analyses presented in the paper.

IMPORTANT: The candidate study itself is protected to due privacy concerns. Interested researchers can however acquire the data themselves from GESIS here https://search.gesis.org/research_data/ZA7704

# Setup

```{r setup, include=FALSE}
# This is the setupchunk, it loads all important packages etc
knitr::opts_chunk$set(echo = TRUE)

# The next bit is quite powerful and useful. 
# First you define which packages you need for your analysis and assign it to 
# the p_needed object. 
p_needed <-
  c("viridis", "stargazer", "MASS", "haven", "ggplot2", "tidyverse", "rmarkdown", "knitr", "lubridate", "janitor", "cowplot", "stringr", "modelsummary", "lmerTest", "plyr", "openxlsx", "mde", "nnet", "viridis")

# Now you check which packages are already installed on your computer.
# The function installed.packages() returns a vector with all the installed 
# packages.
packages <- rownames(installed.packages())
# Then you check which of the packages you need are not installed on your 
# computer yet. Essentially you compare the vector p_needed with the vector
# packages. The result of this comparison is assigned to p_to_install.
p_to_install <- p_needed[!(p_needed %in% packages)]
# If at least one element is in p_to_install you then install those missing
# packages.
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}
# Now that all packages are installed on the computer, you can load them for
# this project. Additionally the expression returns whether the packages were
# successfully loaded.
sapply(p_needed, require, character.only = TRUE)

# This is an option for stargazer tables
# It automatically adapts the output to html or latex,
# depending on whether we want a html or pdf file
stargazer_opt <- ifelse(knitr::is_latex_output(), "latex", "html")
```

# Prep Data

```{r}
can21 <- read_dta("ZA7704_v2-0-0.dta")

can21 <- can21 %>% dplyr::select(d1, # voter vs. party
                                 d2, # own vs. voter
                                 d3, # own vs. party
                          geburtsjahr, # Geburtsjahr
                          partei, #, # Partei
                          geschlecht, # Geschlecht
                          a3f, # prior offices
                          a3g, 
                          a3h,
                          a3i,
                          a3j,
                          a3k,
                          a2a, # prior candidacies
                          a2b, 
                          a2c,
                          a2d,
                          a2e, 
                          gewaehlt # successfull candidacy?
                          )

can21 <- recode_as_na_for(can21, value = 0, criteria = "lt")
```

## New variables

```{r}
# Generate variable who should have priority
can21$rep_prio <- case_when(can21$d1 == 1 & can21$d2 == 2 ~ "voter",
                            can21$d2 == 1 & can21$d3 == 1 ~ "own",
                            can21$d1 == 2 & can21$d3 == 2 ~ "party",
                            can21$d1 == 1 & can21$d2 == 1 & can21$d3 == 2 ~ "circular",
                            can21$d1 == 2 & can21$d2 == 2 & can21$d3 == 1 ~ "circular",
                            .default = NA)

table(can21$rep_prio)

# Age
can21$age <- 2021 - can21$geburtsjahr

# Sub35 age cutoff
can21$sub35 <- ifelse(can21$geburtsjahr > 1986, 1, 0)

# party dummies

can21$p_3 <- ifelse(can21$partei == 3, 1, 0)
can21$p_4 <- ifelse(can21$partei == 4, 1, 0)
can21$p_5 <- ifelse(can21$partei == 5, 1, 0)
can21$p_6 <- ifelse(can21$partei == 6, 1, 0)
can21$p_7 <- ifelse(can21$partei == 7, 1, 0)
can21$p_322 <- ifelse(can21$partei == 322, 1, 0)

# Prior office
can21$prior_office <- ifelse(can21$a3f + can21$a3g + can21$a3h + can21$a3i + can21$a3j + can21$a3k > 0, 1, 0)
table(can21$prior_office)


# Incumbent
can21$incumbent <- ifelse(can21$a2a == 3, 1, 0)
table(can21$incumbent)

# gender to 0/1
can21$geschlecht <- can21$geschlecht - 1
```

# Descriptives

```{r}
can21_sub <- can21 %>% filter(rep_prio %in% c("voter", "own", "party"))

mean(can21_sub$age)

can21_sub %>% filter(prior_office == 0) %>% summarize(age = mean(age))
can21_sub %>% filter(prior_office == 1) %>% summarize(age = mean(age))

can21_sub %>% filter(incumbent == 0) %>% summarize(age = mean(age))
can21_sub %>% filter(incumbent == 1) %>% summarize(age = mean(age))

summary(can21_sub)
```


# Models

```{r}
can21_sub <- can21 %>% filter(rep_prio %in% c("voter", "own", "party"))

model1 <- multinom(rep_prio ~ age + geschlecht + gewaehlt + as.factor(partei),
                   data = can21_sub)
summary(model1)

model2 <- multinom(rep_prio ~ age + incumbent + prior_office + geschlecht + gewaehlt + as.factor(partei),
                   data = can21_sub)
summary(model2)

stargazer(model1, model2,
          covariate.labels = c("Age",
                               "Incumbent",
                               "Prior Office",
                               "Female",
                               "Successful Candidacy",
                               "Party: CSU",
                               "Party: SPD",
                               "Party: FDP",
                               "Party: Greens",
                               "Party: Die LINKE",
                               "Party: AfD",
                               "Intercept"),
          notes = "Baseline for Parties: CDU")
```

# Simulation

## Without Experience

```{r}
# set up
set.seed(06022025)
mu <- coef(model2_ger)
mu <- c(mu[1,], mu[2,])
varcov <- vcov(model2_ger)

nsim <- 1000
J <- 3 # Number of response options
k <- 8 # Number of IVs



X <- as.matrix(cbind(1, can21_sub$age, can21_sub$geschlecht, can21_sub$gewaehlt, can21_sub$p_3, can21_sub$p_4, can21_sub$p_5, can21_sub$p_6, can21_sub$p_7, can21_sub$p_322))

# Sampling distribution
S <- mvrnorm(nsim, mu, varcov)

sel <- 2 # location of age

ageScenario <- seq(min(can21_sub$age), max(can21_sub$age), length.out = 25)

cases <- array(NA, c(dim(X), length(ageScenario)))

cases[, , ] <- X

for (i in 1:length(ageScenario)) {
  cases[, sel, i] <- ageScenario[i]
}

# Calculate the utilities V
V <- array(NA, c(nrow(X), nsim, length(ageScenario), 3))

# Loop over the scenarios
for (i in 1:length(ageScenario)) {
  V[, , i, 1] <-
    apply(matrix(0, nrow = nsim, ncol = ncol(X)), 1, function(s)
      cases[, , i] %*% s)
  V[, , i, 2] <- apply(S[, 1:10], 1, function(s)
    cases[, , i] %*% s)
  V[, , i, 3] <- apply(S[, 11:20], 1, function(s)
    cases[, , i] %*% s)
}

# Summarize over dimensions
Sexp <- apply(V, c(1, 2, 3), function(x)
  sum(exp(x)))

# Get P
P <- array(NA, c(nsim, J, length(ageScenario)))

for (h in 1:length(ageScenario)) {
  for (i in 1:J) {
    P[, i, h] <- apply(exp(V[, , h, i]) / Sexp[, , h], 2, mean)
  }
}

QoImean <- apply(P, c(2, 3), mean)
QoICI <-
  apply(P, c(2, 3), function(x)
    quantile(x, probs = c(0.05, 0.95)))

plot_1 <- ggplot()+
  geom_line(aes(x = ageScenario, y = QoImean[1,], color = "own")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI[1, 1,], ymax = QoICI[2, 1,]), color = plasma(3)[3], fill = plasma(3)[3], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean[2,], color = "party")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI[1, 2,], ymax = QoICI[2, 2,]), color = plasma(3)[2], fill = plasma(3)[2], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean[3,], color = "voter")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI[1, 3,], ymax = QoICI[2, 3,]), color = plasma(3)[1], fill = plasma(3)[1], alpha = 0.3) +
  theme_bw() +
  labs(title = "Excluding Incumbency & Prior Experience",
       x = "Age",
       y = "Predicted Probability",
       caption = "") +
  scale_color_manual(name = "Legend",
                     labels = c("Own Views", "Party's Views", "Voters' Views"),
                     values = c("own" = plasma(3)[3], "party" = plasma(3)[2], "voter" = plasma(3)[1])) +
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 0.9))

# ggsave(file = "pred_prob_1.png", width = 11, height = 7)
```

## With Experience

```{r}
# set up
set.seed(06022025)
mu <- coef(model3_ger)
mu <- c(mu[1,], mu[2,])
varcov <- vcov(model3_ger)

nsim <- 1000
J <- 3 # Number of response options
k <- 10 # Number of IVs


X <- cbind(1, can21_sub$age, can21_sub$incumbent, can21_sub$prior_office, can21_sub$geschlecht, can21_sub$gewaehlt, can21_sub$p_3, can21_sub$p_4, can21_sub$p_5, can21_sub$p_6, can21_sub$p_7, can21_sub$p_322)
X <- X %>% na.omit()
X <- as.matrix(X)

# Sampling distribution
S <- mvrnorm(nsim, mu, varcov)

sel <- 2 # location of age

ageScenario <- seq(min(can21_sub$age), max(can21_sub$age), length.out = 25)

cases <- array(NA, c(dim(X), length(ageScenario)))

cases[, , ] <- X

for (i in 1:length(ageScenario)) {
  cases[, sel, i] <- ageScenario[i]
}

# Calculate the utilities V
V <- array(NA, c(nrow(X), nsim, length(ageScenario), 3))

# Loop over the scenarios
for (i in 1:length(ageScenario)) {
  V[, , i, 1] <-
    apply(matrix(0, nrow = nsim, ncol = ncol(X)), 1, function(s)
      cases[, , i] %*% s)
  V[, , i, 2] <- apply(S[, 1:12], 1, function(s)
    cases[, , i] %*% s)
  V[, , i, 3] <- apply(S[, 13:24], 1, function(s)
    cases[, , i] %*% s)
}

# Summarize over dimensions
Sexp <- apply(V, c(1, 2, 3), function(x)
  sum(exp(x)))

# Get P
P <- array(NA, c(nsim, J, length(ageScenario)))

for (h in 1:length(ageScenario)) {
  for (i in 1:J) {
    P[, i, h] <- apply(exp(V[, , h, i]) / Sexp[, , h], 2, mean)
  }
}

QoImean_2 <- apply(P, c(2, 3), mean)
QoICI_2 <-
  apply(P, c(2, 3), function(x)
    quantile(x, probs = c(0.05, 0.95)))


plot_2 <- ggplot()+
  geom_line(aes(x = ageScenario, y = QoImean_2[1,], color = "own")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 1,], ymax = QoICI_2[2, 1,]), color = plasma(3)[3], fill = plasma(3)[3], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean_2[2,], color = "party")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 2,], ymax = QoICI_2[2, 2,]), color = plasma(3)[2], fill = plasma(3)[2], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean_2[3,], color = "voter")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 3,], ymax = QoICI_2[2, 3,]), color = plasma(3)[1], fill = plasma(3)[1], alpha = 0.3) +
  theme_bw() +
  labs(title = "Including Incumbency & Prior Experience",
       x = "Age",
       y = "",
       caption = "Shaded Areas = 90% CIs") +
  scale_color_manual(name = "Legend",
                     labels = c("Own Views", "Party's Views", "Voters' Views"),
                     values = c("own" = plasma(3)[3], "party" = plasma(3)[2], "voter" = plasma(3)[1])) +
  scale_y_continuous(limits = c(0, 0.9))

legend <- get_legend(plot_2)

plot_2 <- ggplot()+
  geom_line(aes(x = ageScenario, y = QoImean_2[1,], color = "own")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 1,], ymax = QoICI_2[2, 1,]), color = plasma(3)[3], fill = plasma(3)[3], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean_2[2,], color = "party")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 2,], ymax = QoICI_2[2, 2,]), color = plasma(3)[2], fill = plasma(3)[2], alpha = 0.3) +
  geom_line(aes(x = ageScenario, y = QoImean_2[3,], color = "voter")) +
  geom_ribbon(aes(x = ageScenario, ymin = QoICI_2[1, 3,], ymax = QoICI_2[2, 3,]), color = plasma(3)[1], fill = plasma(3)[1], alpha = 0.3) +
  theme_bw() +
  labs(title = "Including Incumbency & Prior Experience",
       x = "Age",
       y = "",
       caption = "Shaded Areas = 90% CIs") +
  scale_color_manual(name = "Legend",
                     labels = c("Own Views", "Party's Views", "Voters' Views"),
                     values = c("own" = plasma(3)[3], "party" = plasma(3)[2], "voter" = plasma(3)[1])) +
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 0.9))

#ggsave(file = "pred_prob_2.png", width = 11, height = 7)
```

```{r}
grid_1 <- plot_grid(plot_1, plot_2, rel_widths = c(1, 1))
plot_grid(grid_1, legend, ncol = 2, rel_widths = c(1, 0.2))
#ggsave2(file = "plot_grid_RR.png", width = 11, height = 6)
```

# Models only Successfull

```{r}
can21_gew <- can21_sub %>% filter(gewaehlt == 1)

model1_gew <- multinom(rep_prio ~ age + geschlecht + as.factor(partei),
                   data = can21_gew)
summary(model1_gew)

model2_gew <- multinom(rep_prio ~ age + incumbent + prior_office + geschlecht + as.factor(partei),
                   data = can21_gew)
summary(model2_gew)

stargazer(model1_ger, model2_ger,
          covariate.labels = c("Age",
                               "Incumbent",
                               "Prior Office",
                               "Female",
                               "Successful Candidacy",
                               "Party: CSU",
                               "Party: SPD",
                               "Party: FDP",
                               "Party: Greens",
                               "Party: Die LINKE",
                               "Party: AfD",
                               "Intercept"),
          notes = "Baseline for Parties: CDU")
```